home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / ztrsm.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  13.6 KB  |  418 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  5.      $                   B, LDB )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
  8.       INTEGER            M, N, LDA, LDB
  9.       COMPLEX*16         ALPHA
  10. *     .. Array Arguments ..
  11.       COMPLEX*16         A( LDA, * ), B( LDB, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  ZTRSM  solves one of the matrix equations
  18. *
  19. *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
  20. *
  21. *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  22. *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  23. *
  24. *     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ).
  25. *
  26. *  The matrix X is overwritten on B.
  27. *
  28. *  Parameters
  29. *  ==========
  30. *
  31. *  SIDE   - CHARACTER*1.
  32. *           On entry, SIDE specifies whether op( A ) appears on the left
  33. *           or right of X as follows:
  34. *
  35. *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
  36. *
  37. *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
  38. *
  39. *           Unchanged on exit.
  40. *
  41. *  UPLO   - CHARACTER*1.
  42. *           On entry, UPLO specifies whether the matrix A is an upper or
  43. *           lower triangular matrix as follows:
  44. *
  45. *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  46. *
  47. *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  48. *
  49. *           Unchanged on exit.
  50. *
  51. *  TRANSA - CHARACTER*1.
  52. *           On entry, TRANSA specifies the form of op( A ) to be used in
  53. *           the matrix multiplication as follows:
  54. *
  55. *              TRANSA = 'N' or 'n'   op( A ) = A.
  56. *
  57. *              TRANSA = 'T' or 't'   op( A ) = A'.
  58. *
  59. *              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
  60. *
  61. *           Unchanged on exit.
  62. *
  63. *  DIAG   - CHARACTER*1.
  64. *           On entry, DIAG specifies whether or not A is unit triangular
  65. *           as follows:
  66. *
  67. *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  68. *
  69. *              DIAG = 'N' or 'n'   A is not assumed to be unit
  70. *                                  triangular.
  71. *
  72. *           Unchanged on exit.
  73. *
  74. *  M      - INTEGER.
  75. *           On entry, M specifies the number of rows of B. M must be at
  76. *           least zero.
  77. *           Unchanged on exit.
  78. *
  79. *  N      - INTEGER.
  80. *           On entry, N specifies the number of columns of B.  N must be
  81. *           at least zero.
  82. *           Unchanged on exit.
  83. *
  84. *  ALPHA  - COMPLEX*16      .
  85. *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
  86. *           zero then  A is not referenced and  B need not be set before
  87. *           entry.
  88. *           Unchanged on exit.
  89. *
  90. *  A      - COMPLEX*16       array of DIMENSION ( LDA, k ), where k is m
  91. *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
  92. *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
  93. *           upper triangular part of the array  A must contain the upper
  94. *           triangular matrix  and the strictly lower triangular part of
  95. *           A is not referenced.
  96. *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
  97. *           lower triangular part of the array  A must contain the lower
  98. *           triangular matrix  and the strictly upper triangular part of
  99. *           A is not referenced.
  100. *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
  101. *           A  are not referenced either,  but are assumed to be  unity.
  102. *           Unchanged on exit.
  103. *
  104. *  LDA    - INTEGER.
  105. *           On entry, LDA specifies the first dimension of A as declared
  106. *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
  107. *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
  108. *           then LDA must be at least max( 1, n ).
  109. *           Unchanged on exit.
  110. *
  111. *  B      - COMPLEX*16       array of DIMENSION ( LDB, n ).
  112. *           Before entry,  the leading  m by n part of the array  B must
  113. *           contain  the  right-hand  side  matrix  B,  and  on exit  is
  114. *           overwritten by the solution matrix  X.
  115. *
  116. *  LDB    - INTEGER.
  117. *           On entry, LDB specifies the first dimension of B as declared
  118. *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  119. *           max( 1, m ).
  120. *           Unchanged on exit.
  121. *
  122. *
  123. *  Level 3 Blas routine.
  124. *
  125. *  -- Written on 8-February-1989.
  126. *     Jack Dongarra, Argonne National Laboratory.
  127. *     Iain Duff, AERE Harwell.
  128. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  129. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  130. *
  131. *
  132. *     .. External Functions ..
  133.       LOGICAL            LSAME
  134.       EXTERNAL           LSAME
  135. *     .. External Subroutines ..
  136.       EXTERNAL           XERBLA
  137. *     .. Intrinsic Functions ..
  138.       INTRINSIC          DCONJG, MAX
  139. *     .. Local Scalars ..
  140.       LOGICAL            LSIDE, NOCONJ, NOUNIT, UPPER
  141.       INTEGER            I, INFO, J, K, NROWA
  142.       COMPLEX*16         TEMP
  143. *     .. Parameters ..
  144.       COMPLEX*16         ONE
  145.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  146.       COMPLEX*16         ZERO
  147.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  148. *     ..
  149. *     .. Executable Statements ..
  150. *
  151. *     Test the input parameters.
  152. *
  153.       LSIDE  = LSAME( SIDE  , 'L' )
  154.       IF( LSIDE )THEN
  155.          NROWA = M
  156.       ELSE
  157.          NROWA = N
  158.       END IF
  159.       NOCONJ = LSAME( TRANSA, 'T' )
  160.       NOUNIT = LSAME( DIAG  , 'N' )
  161.       UPPER  = LSAME( UPLO  , 'U' )
  162. *
  163.       INFO   = 0
  164.       IF(      ( .NOT.LSIDE                ).AND.
  165.      $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
  166.          INFO = 1
  167.       ELSE IF( ( .NOT.UPPER                ).AND.
  168.      $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
  169.          INFO = 2
  170.       ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
  171.      $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
  172.      $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
  173.          INFO = 3
  174.       ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
  175.      $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
  176.          INFO = 4
  177.       ELSE IF( M  .LT.0               )THEN
  178.          INFO = 5
  179.       ELSE IF( N  .LT.0               )THEN
  180.          INFO = 6
  181.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  182.          INFO = 9
  183.       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
  184.          INFO = 11
  185.       END IF
  186.       IF( INFO.NE.0 )THEN
  187.          CALL XERBLA( 'ZTRSM ', INFO )
  188.          RETURN
  189.       END IF
  190. *
  191. *     Quick return if possible.
  192. *
  193.       IF( N.EQ.0 )
  194.      $   RETURN
  195. *
  196. *     And when  alpha.eq.zero.
  197. *
  198.       IF( ALPHA.EQ.ZERO )THEN
  199.          DO 20, J = 1, N
  200.             DO 10, I = 1, M
  201.                B( I, J ) = ZERO
  202.    10       CONTINUE
  203.    20    CONTINUE
  204.          RETURN
  205.       END IF
  206. *
  207. *     Start the operations.
  208. *
  209.       IF( LSIDE )THEN
  210.          IF( LSAME( TRANSA, 'N' ) )THEN
  211. *
  212. *           Form  B := alpha*inv( A )*B.
  213. *
  214.             IF( UPPER )THEN
  215.                DO 60, J = 1, N
  216.                   IF( ALPHA.NE.ONE )THEN
  217.                      DO 30, I = 1, M
  218.                         B( I, J ) = ALPHA*B( I, J )
  219.    30                CONTINUE
  220.                   END IF
  221.                   DO 50, K = M, 1, -1
  222.                      IF( B( K, J ).NE.ZERO )THEN
  223.                         IF( NOUNIT )
  224.      $                     B( K, J ) = B( K, J )/A( K, K )
  225.                         DO 40, I = 1, K - 1
  226.                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
  227.    40                   CONTINUE
  228.                      END IF
  229.    50             CONTINUE
  230.    60          CONTINUE
  231.             ELSE
  232.                DO 100, J = 1, N
  233.                   IF( ALPHA.NE.ONE )THEN
  234.                      DO 70, I = 1, M
  235.                         B( I, J ) = ALPHA*B( I, J )
  236.    70                CONTINUE
  237.                   END IF
  238.                   DO 90 K = 1, M
  239.                      IF( B( K, J ).NE.ZERO )THEN
  240.                         IF( NOUNIT )
  241.      $                     B( K, J ) = B( K, J )/A( K, K )
  242.                         DO 80, I = K + 1, M
  243.                            B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
  244.    80                   CONTINUE
  245.                      END IF
  246.    90             CONTINUE
  247.   100          CONTINUE
  248.             END IF
  249.          ELSE
  250. *
  251. *           Form  B := alpha*inv( A' )*B
  252. *           or    B := alpha*inv( conjg( A' ) )*B.
  253. *
  254.             IF( UPPER )THEN
  255.                DO 140, J = 1, N
  256.                   DO 130, I = 1, M
  257.                      TEMP = ALPHA*B( I, J )
  258.                      IF( NOCONJ )THEN
  259.                         DO 110, K = 1, I - 1
  260.                            TEMP = TEMP - A( K, I )*B( K, J )
  261.   110                   CONTINUE
  262.                         IF( NOUNIT )
  263.      $                     TEMP = TEMP/A( I, I )
  264.                      ELSE
  265.                         DO 120, K = 1, I - 1
  266.                            TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
  267.   120                   CONTINUE
  268.                         IF( NOUNIT )
  269.      $                     TEMP = TEMP/DCONJG( A( I, I ) )
  270.                      END IF
  271.                      B( I, J ) = TEMP
  272.   130             CONTINUE
  273.   140          CONTINUE
  274.             ELSE
  275.                DO 180, J = 1, N
  276.                   DO 170, I = M, 1, -1
  277.                      TEMP = ALPHA*B( I, J )
  278.                      IF( NOCONJ )THEN
  279.                         DO 150, K = I + 1, M
  280.                            TEMP = TEMP - A( K, I )*B( K, J )
  281.   150                   CONTINUE
  282.                         IF( NOUNIT )
  283.      $                     TEMP = TEMP/A( I, I )
  284.                      ELSE
  285.                         DO 160, K = I + 1, M
  286.                            TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J )
  287.   160                   CONTINUE
  288.                         IF( NOUNIT )
  289.      $                     TEMP = TEMP/DCONJG( A( I, I ) )
  290.                      END IF
  291.                      B( I, J ) = TEMP
  292.   170             CONTINUE
  293.   180          CONTINUE
  294.             END IF
  295.          END IF
  296.       ELSE
  297.          IF( LSAME( TRANSA, 'N' ) )THEN
  298. *
  299. *           Form  B := alpha*B*inv( A ).
  300. *
  301.             IF( UPPER )THEN
  302.                DO 230, J = 1, N
  303.                   IF( ALPHA.NE.ONE )THEN
  304.                      DO 190, I = 1, M
  305.                         B( I, J ) = ALPHA*B( I, J )
  306.   190                CONTINUE
  307.                   END IF
  308.                   DO 210, K = 1, J - 1
  309.                      IF( A( K, J ).NE.ZERO )THEN
  310.                         DO 200, I = 1, M
  311.                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  312.   200                   CONTINUE
  313.                      END IF
  314.   210             CONTINUE
  315.                   IF( NOUNIT )THEN
  316.                      TEMP = ONE/A( J, J )
  317.                      DO 220, I = 1, M
  318.                         B( I, J ) = TEMP*B( I, J )
  319.   220                CONTINUE
  320.                   END IF
  321.   230          CONTINUE
  322.             ELSE
  323.                DO 280, J = N, 1, -1
  324.                   IF( ALPHA.NE.ONE )THEN
  325.                      DO 240, I = 1, M
  326.                         B( I, J ) = ALPHA*B( I, J )
  327.   240                CONTINUE
  328.                   END IF
  329.                   DO 260, K = J + 1, N
  330.                      IF( A( K, J ).NE.ZERO )THEN
  331.                         DO 250, I = 1, M
  332.                            B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
  333.   250                   CONTINUE
  334.                      END IF
  335.   260             CONTINUE
  336.                   IF( NOUNIT )THEN
  337.                      TEMP = ONE/A( J, J )
  338.                      DO 270, I = 1, M
  339.                        B( I, J ) = TEMP*B( I, J )
  340.   270                CONTINUE
  341.                   END IF
  342.   280          CONTINUE
  343.             END IF
  344.          ELSE
  345. *
  346. *           Form  B := alpha*B*inv( A' )
  347. *           or    B := alpha*B*inv( conjg( A' ) ).
  348. *
  349.             IF( UPPER )THEN
  350.                DO 330, K = N, 1, -1
  351.                   IF( NOUNIT )THEN
  352.                      IF( NOCONJ )THEN
  353.                         TEMP = ONE/A( K, K )
  354.                      ELSE
  355.                         TEMP = ONE/DCONJG( A( K, K ) )
  356.                      END IF
  357.                      DO 290, I = 1, M
  358.                         B( I, K ) = TEMP*B( I, K )
  359.   290                CONTINUE
  360.                   END IF
  361.                   DO 310, J = 1, K - 1
  362.                      IF( A( J, K ).NE.ZERO )THEN
  363.                         IF( NOCONJ )THEN
  364.                            TEMP = A( J, K )
  365.                         ELSE
  366.                            TEMP = DCONJG( A( J, K ) )
  367.                         END IF
  368.                         DO 300, I = 1, M
  369.                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
  370.   300                   CONTINUE
  371.                      END IF
  372.   310             CONTINUE
  373.                   IF( ALPHA.NE.ONE )THEN
  374.                      DO 320, I = 1, M
  375.                         B( I, K ) = ALPHA*B( I, K )
  376.   320                CONTINUE
  377.                   END IF
  378.   330          CONTINUE
  379.             ELSE
  380.                DO 380, K = 1, N
  381.                   IF( NOUNIT )THEN
  382.                      IF( NOCONJ )THEN
  383.                         TEMP = ONE/A( K, K )
  384.                      ELSE
  385.                         TEMP = ONE/DCONJG( A( K, K ) )
  386.                      END IF
  387.                      DO 340, I = 1, M
  388.                         B( I, K ) = TEMP*B( I, K )
  389.   340                CONTINUE
  390.                   END IF
  391.                   DO 360, J = K + 1, N
  392.                      IF( A( J, K ).NE.ZERO )THEN
  393.                         IF( NOCONJ )THEN
  394.                            TEMP = A( J, K )
  395.                         ELSE
  396.                            TEMP = DCONJG( A( J, K ) )
  397.                         END IF
  398.                         DO 350, I = 1, M
  399.                            B( I, J ) = B( I, J ) - TEMP*B( I, K )
  400.   350                   CONTINUE
  401.                      END IF
  402.   360             CONTINUE
  403.                   IF( ALPHA.NE.ONE )THEN
  404.                      DO 370, I = 1, M
  405.                         B( I, K ) = ALPHA*B( I, K )
  406.   370                CONTINUE
  407.                   END IF
  408.   380          CONTINUE
  409.             END IF
  410.          END IF
  411.       END IF
  412. *
  413.       RETURN
  414. *
  415. *     End of ZTRSM .
  416. *
  417.       END
  418.